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Abstract 

With the growing use of computer modeling 
and simulation, in all aspects of engineering, the 
scope of traditional optimization has to be ex- 
tended to include simulation models. Some unique 
aspects have to be addressed while optimizing via 
stochastic simulation models. The optimization 
procedure has to explicitly account for the ran- 
domness inherent in the stochastic measures pre- 
dicted by the model. This paper outlines a general- 
purpose framework for optimization of terminat- 
ing discrete-event simulation models. The meth- 
odology combines a chance constraint approach 
for problem formulation, together with standard 
statistical estimation and analyses techniques. The 
applicability of the optimization framework is il- 
lustrated by minimizing the operation and support 
resources of a launch vehicle, through a simulation 
model. 

Introduction 

Simulations are used to study physical systems 
or processes that are too complex to permit ana- 
lytical model formulation and evaluation. The 
complexity in the physical system is usually at- 
tributed to its size, interaction between its sub- 
systems or components, and randomness inherent 
in its processes. In simulation modeling, the 
physical system under consideration is mimicked 


or imitated with the help of a computer program. 
The performance or the response of the system 
can be observed by running the simulation pro- 
gram with operational data. These operating pa- 
rameters can be changed to obtain a desired im- 
provement in the simulation response and hence 
performance of the system. 

The power of modern computing and the 
emergence of simulation software in recent years, 
has contributed to the increasing use of simulation 
in industry. Simulation models now span diverse 
application areas such as: designing and analyzing 
manufacturing systems, characterization of semi- 
conductor devices, determining ordering policies 
for an inventory system, and analyzing a financial 
or economic system, to name a few (Law and 
Kelton, 1991). 

With the growing incidence of simulation mod- 
eling, it is essential to extend the scope of tradi- 
tional optimization to include the simulation do- 
main. Optimization procedures for simulation 
have to specifically account for the randomness 
inherent in the stochastic simulation models. Often 
in practice, however, the randomness is ignored 
and simulation models are analyzed in an ad-hoc 
manner (Law and Kelton, 1991). In this paper, a 
unifying framework, based on well established 
procedures in mathematical programming and sta- 
tistics, is outlined for use with terminating types of 
discrete-event simulations. Such a methodology 
has been generally lacking in the simulation opti- 



mization literature, which largely focuses on de- 
veloping new optimization approaches. 

The framework employs the chance constraint 
approach for treating stochastic constraints for the 
purpose of problem formulation (Charnes and 
Cooper, 1956). Standard statistical techniques are 
employed for estimating the stochastic compo- 
nents of the simulation. Optimization techniques 
such as path or pattern search methods, can then 
be applied in conjunction with this framework, to 
optimize the system. The framework provides 
guidelines for a rigorous study so that resulting 
inferences can be made with statistical confidence. 
The methodology is illustrated by minimizing the 
operation and support resources for a reusable 
launch vehicle, modeled through a terminating dis- 
crete-event simulation. A genetic algorithm is used 
to perform the numerical optimization. 

The Simulation Optimization Problem 

In this section, the general stochastic simulation 
optimization problem is introduced with the help 
of an example. Consider a terminating simulation 
modeling the activities of a bank telling system, in 
an optimization study. The model simulates the 
essentially stochastic processes of customer arri- 
val, queues, and teller service and idle times, 
through appropriate probability distributions. Let 
us suppose that the optimization problem is as 
follows: 

Minimize the average waiting time for bank 
customers arriving randomly, while ensuring 
that the average teller idle time does not exceed a 
certain level (say 20 minutes), and resources re- 
quired do not exceed an allowable limit (say 
$60,000/ year). 

The objective function i.e. ‘minimization of the 
average time that a customer waits for service’, is 
stochastic in nature and predicted by the simula- 
tion. The problem has two constraints. The 
‘average teller idle time not to exceed 20 minutes’ 
constraint is stochastic and predicted by the simu- 


lation. The ‘total resources (teller salaries) not to 
exceed $60,000’ constraint is deterministic. Sala- 
ries are based on the number of tellers assigned to 
the simulation by the user at run-time. The tellers 
assigned clearly influence the resources and sys- 
tem performance. Optimization thus requires de- 
termining the number of tellers to assign, so as to 
minimize the average customer waiting time 
within the constraint space. 

A general simulation optimization problem can 
thus be stated as determining the settings of the 

input parameters X = (Xj, Xj„ , Xp,) e S, 

over a region S cz RP , which optimize the ex- 
pected response E(F(X)) (Jacobson and Schruben, 
1989). The region S cr RP is defined or bounded 
by the constraints of the form E(Z(X)) > b, where 
Z(X) may be based on the stochastic simulation 
response. 

Such problems are known to be hard to solve 
(Jacobson and Schruben, 1989). The objective 
function and/or constraints may be predicted by 
the stochastic simulation and hence the problem 
lacks a closed form analytical formulation. Calcu- 
lus based methods and mathematical program- 
ming techniques, requiring explicit problem for- 
mulation and continuous domains, cannot be ap- 
plied (Fu, 1994). Non-linear programming ap- 
proaches such as path and pattern search methods 
are generally used to optimize simulated systems. 
Recently, probabilistic search strategies such as 
genetic algorithms and simulated annealing have 
also been used with good results. 

Assumptions 

Due to its random nature, each run of a simula- 
tion has the potential to give rise to slightly differ- 
ent output responses. Hence, the responses pre- 
dicted by the simulation are generally estimated 
through a sample collected by replicating the 
simulation with different random number seeds. 
Statistics such as the mean, variance, and interval 
estimates are used for estimation. These statistics, 
computed in the standard manner, are based on 
several assumptions. These are: 



1. The stochastic process is covariance stationary. 

2. The sample variance is an unbiased estimator of 

the population variance. 

3. The observations are independent and identi- 
cally distributed. 

The above assumptions do not always hold true in 
simulation studies. For example, covariance sta- 
tionarity may not rigorously hold for terminating 
simulations, unless the simulation time span is 
sufficiently long to warrant stationarity. In termi- 
nating simulations, a simulation stops when a 
natural events signaling the end of the simulation 
run occurs. (For instance, in the bank teller exam- 
ple, the simulation time span may be defined as 8 
a.m. to 5 p.m.). Since we do not necessarily run 
the simulation until steady state, the underlying 
joint distributions of the random variables may 
change over time. Also, in order to generate inde- 
pendent and identically distributed observations, a 
true random number generator is required. How- 
ever, in practice, pseudo-random number genera- 
tors are used. Furthermore, Law and Kelton 
(1991) observe that simulation output data is usu- 
ally correlated. 

Therefore, the above assumptions do not hold 
in the strictest sense and may be violated by a 
varying degree. However, there are no other ac- 
cepted analyses methods for simulation data. 
Hence it is recommended that regardless of slight 
violations, standard statistical estimation be used 
(Law and Kelton, 1991; Pritsker, 1984; Kleijnen, 
1987; and Fishman, 1978). 

Accuracy 

One of the issues associated with a stochastic 
simulation, is the accuracy with which a stochastic 
variable is to be estimated. The desired accuracy 
for each parameter can be specified by the decision 
maker in terms of statistical confidence intervals. 
(Confidence intervals state the probability ( l-a), 
that the true mean is actually contained in an inter- 
val of width (±w), about the estimated mean). 

When a simulation involves multiple stochastic 
variables, the overall confidence (l-a ) associated 


with an optimization study, is based on satisfying 
the individual confidence intervals (l-a [) simul- 
taneously. Thus the overall confidence satisfies the 
Bonferroni inequality given by 

p> l-La,, 

and implies a lower overall accuracy. 

For studies involving ten or less stochastic vari- 
ables, if an overall confidence for the study (l-a ) 
is desired, then the individual confidences (l-a , ) 
can be selected by the following relation: 

L«i = «- 

However, for more than ten stochastic variables, 
the accuracy of individual variables obtained by the 
Bonferroni inequality may be prohibitively high. 
For example, if an overall confidence of 90% is 
desired for a simulation involving ten variables, 
the individual confidences have to be atleast 99%. 
(Similarly, if the individual confidence intervals of 
ten stochastic variables is 90%, the overall confi- 
dence is only greater than or equal to zero). There- 
fore, for such cases, standard 90% or 95% individ- 
ual confidence intervals are recommended (Law 
and Kelton, 1991). However, the analysis results 
in such cases should be interpreted with caution, as 
one or more of the individual confidence intervals 
may not contain the corresponding true mean. 

Replications 

Once the accuracy for each stochastic variable is 
established, the number of replications required for 
the optimization study can be determined. The ac- 
curacy depends on the sample size or the number 
of replications used in a study. A large sample size 
implies that the estimated mean is closer to the true 
mean, and hence increases the accuracy. A high 
level of accuracy is usually desired so that the re- 
sults of the simulation study and hence subsequent 
decisions can be made with a satisfactory level of 
confidence. However, due to the finite resources 
(CPU time) available for a simulation study, the 



number of replications possible are usually lim- 
ited. 

The individual confidence levels established 
above can be used to obtain the required sample 
size for a desired accuracy in estimating the indi- 
vidual means. The following steps are undertaken 
to compute the sample size (Law and Kelton, 
1991; Kleijnen, 1987): 

1. In an initial pilot experiment, the simulation 
model is run for say 1000 or more replications, 
to obtain a representative distribution for each 
stochastic variable. 

2. The estimated mean and variance for each sto- 
chastic parameter are computed from the pilot 
experiment. 

3. Based on these estimates, the sample size or the 
required number of replications can be com- 
puted for each stochastic variable. For the speci- 
fied confidence or probability (1 -a J that an 
interval (width ±w about the estimated mean) 
contains the true mean, the number of replica- 
tions is determined as follows: 

fe-B S / W X 

where z,__ a represents the standard normal sta- 
tistic covering an area (1 -a ), and 5 is the stan- 
dard deviation of the sample. 

The highest number of replications so computed is 
selected for the simulation study. The individual 
confidence of variables can be assessed at this new 
sample size. The variables that required lower rep- 
lications, will have correspondingly higher confi- 
dences at the newly established sample size. 

Chance Constraints 

The variability inherent in stochastic constraints 
complicates the simulation optimization problem 
by forming fuzzy boundaries for the feasible re- 
gion. This presents a danger of erroneously ac- 


cepting a solution as feasible, while it may have a 
high probability of being infeasible, and vice versa. 
The chance constraint theory approach can be used 
to convert the stochastic constraints into determi- 
nistic constraints. Chance constraints (Chames and 
Cooper, 1959) permit constraint violation up to a 
pre-specified probability limit. The decision maker 
expresses a risk tolerance, in terms of a permissi- 
ble probability of constraint violation. Consider a 
stochastic constraint of the form A(x) < b, where 
A(x) is a simulation response. Using the chance 
constraint approach, this can be reformulated in 
terms of risk tolerance as P(A(x) > b) < a, where 
a denotes the extent to which constraint violation 
is permitted. 

The chance constraints can be implemented 
through confidence interval estimates (Teleb and 
Azadivar, 1994). We know that the confidence 
associated with an interval estimate denotes the 
probability that the interval contains the true pa- 
rameter. For example, the upper limit associated 
with an interval of confidence (1 -a) states that the 
probability that the true mean exceeds this limit is 
at most a. . Thus, the upper and lower limits of the 
interval at the specified confidence (1 -a) (or risk 
a) provide deterministic boundaries for the in- 
feasible region. In this manner, confidence inter- 
vals provide a means of implementing chance 
constraints. 

Using confidence intervals, the upper limit at 
confidence (1 -a) can be used to denote the chance 
constraint P(A(x) > b) < a. as A(x) upplim a < b. The 
confidence intervals are estimated by using the 
Student’s t distribution in the standard manner: 


— r 

v n 

_ s) 

Low_ Lim= x 

Vn 

where n is the number of replications, x is the es- 
timated mean, s is the standard deviation, and t (n _ 
1)x is obtained from the Student 1 distribution at (n- 
1 ) degrees of freedom and a coverage. 



The above confidence intervals based on Stu- 
dent’s t are robust to minor deviations from nor- 
mality of the distribution. However, in cases of 
serious non-normality andvery small sample size, 
the Johnson’s modified t statistic for non-normal 
distributions, is recommended (Johnson, 1978; 
Kleijnen, 1987). This adjusted statistic approxi- 
mates Student's t distribution, by accounting for 
the skewness of the distribution, thus permitting 
its use for hypothesis testing and confidence inter- 
vals. Johnson's modified t statistic has been used 
successfully on distributions with varying degrees 
of non-normality, including the exponential distri- 
bution (Johnson, 1978; Kleijnen, 1985). The con- 
fidence intervals by the modified statistic are given 
by 


Upp_Lim= x+ = +■ J 
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Low_Lim= x— 


ft-(n-l)/* S ^_|_ 7 3 

t/h 6s 2 n 


where 3 is the third central moment estimated in 
the standard manner. 

Application to Operations and 
Support Modeling 

The above simulation optimization framework 
is applied to a NASA model simulating the opera- 
tion and support activities of launch vehicles dur- 
ing conceptual design (Ebeling and Donohue, 
1994). It uses estimated values for component re- 
liability and maintainability, to simulate the mis- 
sions, and pre- and post-flight maintenance. Given 
the operation and support resources (individual 
crew assigned to nine maintenance sub-systems 
c n , and vehicles in the fleet v) the model predicts 
the number of missions flown and the mean vehi- 
cle tum-time for a particular space program. The 
problem is to: 


subject to stochastic constraints evaluated by the 
simulation: 

(1) Average delay E( delay) not to exceed 48 

hours 

(2) Average mission rate of the space program 
E(suc_mis) is met (140 missions in five 
years). 

Minimization of the resources thus requires the 
determination of the smallest fleet and mainte- 
nance crew size, that enables meeting the target 
mission rate in a timely manner. 

Problem Formulation The stochastic constraints 
were converted to deterministic constraints by us- 
ing the chance constraint approach. The NASA 
decision makers expressed a risk tolerance, in 
terms of permissible probability of constraint vio- 
lation, of 5%. These risk tolerances are stated 
symbolically as: 

P[E( delay) > 48 hours] < 0.05 
P[E(suc_ mis) < 140 ] < 0.05 

The limits at the specified 95% confidence (5% 
risk tolerance) were established by using the Stu- 
dent’s t distribution as described above. These 
limits provide deterministic boundaries for the in- 
feasible region, as follows: 

delay U p p Ji my ,()5 ^ 48 hours 
suc_mis low lim 05 > 140 

Accuracy and Sample Size The accuracy desired 
was specified in terms of confidence (1 -a ) that 
the true mean is within the interval ±w, by NASA 
engineers as: 



Desired 
width ±w 

Desired 
conf. (1-a ) 

Delay 

±48 hours 

80% 

Successful Missions 

±2 

95% 


Minimize the cost function The number of replications required were com- 

100*v+ c j+c 2 +cj+C4+ cj+ctf + cy+cg+ eg puted using the steps outlined in the preceding 



sections. This computation yielded: 19 replications 
for the average delay, and 18 replications for suc- 
cessful missions. Hence, a conservative sample 
size of 20 was selected. 

Results The problem formulated above was 
solved by applying a genetic algorithm (Gage, 
1995). A genetic algorithm is a biased random 
search approach that has been used successfully 
for optimization in complex landscapes (Holland 
1975; Goldberg, 1989). Unlike other mathematical 
programming techniques, a genetic algorithm does 
not require continuity, differentiability, and explicit 
problem formulation. Due to these reasons, the 
genetic algorithm is an attractive technique for the 
simulation domain. 

Inspired by the mechanics of natural selection 
and evolution, a genetic algorithm uses binary 
strings to encode parameters, for representing so- 
lutions in the search space. Based on the fitness of 
such binary strings, contained in a randomly se- 
lected initial generation, subsequent generations are 
created by mating and mutation of the strings. The 
genetic algorithm continues to form new genera- 
tions, in this manner, until a suitable criterion is 
reached. For a detailed description of genetic algo- 
rithms please see Goldberg (1989), and Srinivas 
and Patnaik (1994). 

The minimum objective function found by the 
genetic algorithm when applied to the problem 
formulated above was 273. Two vehicles and 73 
crew contribute to this objective function. The con- 
straints at this input combination were satisfied 
and well within their tolerance levels. The mean of 
average delay was 15.73 hours, with a 95% confi- 
dence that the delay does not exceed 38 hours. The 
target mission rate of 140 missions in a five year 
time span was achieved with a 95% confidence. 
The above stochastic measures were estimated at 
an 80% confidence of being within ± 2 days, and 
95% confidence of being within ± 2 missions, re- 
spectively. The overall accuracy of the optimiza- 
tion study was lower at 75% confidence. 


Conclusions 

Stochastic simulation modeling presents some 
unique problems that have to be addressed during 
optimization. In this paper a general purpose 
framework has been outlined for the optimization 
of terminating discrete event simulations. A uni- 
fying framework or methodology for treating the 
stochastic components of such problems has been 
generally lacking in the literature. The methodol- 
ogy presented here employs a chance constraint 
approach to convert the fuzzy boundaries pre- 
sented by stochastic constraints into deterministic 
constraints. Standard statistical techniques are used 
to perform the optimization study at desired levels 
of confidence. Therefore the optimization frame- 
work is based on a sound statistical background. 
Inferences drawn from such studies can be stated 
in terms of interval estimates, which are more ro- 
bust than point estimates. This provides more in- 
formation while making decisions based on the 
simulation study. The practical significance of the 
methodology was illustrated by applying it to a 
NASA simulation. The operation and support re- 
sources of reusable launch vehicles were success- 
fully minimized with the help of a genetic algo- 
rithm. 
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